#include "FEMSpace_unstruct.h"
#include "function2.h"
#include <ctime>
#include <math.h>
#include <typeinfo>
int main(int argc, char *argv[]) {
    // g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    vector<Point<2>> p = {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}};
    Domain<2> *domain = new Rectangle2D(p);
    int order = atoi(argv[1]);
    EasyMesh *m = new EasyMesh("D", order);
    assert(order >= 1 && order <= 4);
    TemplateElement<2> *e;
    if (order == 1) {
        e = new Triangular_1_Element();
    } else if (order == 2) {
        e = new Triangular_2_Element();
    } else if (order == 3) {
        e = new Triangular_3_Element();
    } else if (order == 4) {
        e = new Triangular_4_Element();
    }
    Equation<2> *equ = new PossionEquation<2>();
    equ->SetRightHandsTermFunction(fe);
    equ->SetBoundaryConditionFunction({g1, g2, g3, g4});
    Possion_2D_Unstruct possionproblem(m, e, equ, domain);
    clock_t start, end;
    start = clock();
    possionproblem.setup_system();
    end = clock(); // 程序结束用时
    double endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "setup_system takes  " << endtime * 1000 << " ms"
         << endl; // ms为单位
    start = clock();
    possionproblem.assemble_system();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "assemble_system takes  " << endtime * 1000 << " ms"
         << endl; // ms为单位
    start = clock();
    possionproblem.dealwithBoundaryCondition();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "dealwithBoundaryCondition takes  " << endtime * 1000 << " ms"
         << endl; // ms为单位
    start = clock();
    possionproblem.solve();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "solve takes  " << endtime * 1000 << " ms" << endl; // ms为单位
    int numdofs = m->getTotalNumDofs();
    Eigen::VectorXd Truesolution = Eigen::MatrixXd::Zero(numdofs, 1);
    vector<Dofs<2>> dofsList = m->getAllDofs();
    for (auto &dof : dofsList) {
        Point<2> pnt = dof.getPoint();
        double point[2] = {pnt[0], pnt[1]};
        int index = dof.getGlobalIndex();
        Truesolution[index] = u(point);
    }

    Eigen::VectorXd residual =
        possionproblem.A() * possionproblem.getSolution() -
        possionproblem.Rhs();
    double residual_l2 = residual.squaredNorm();
    double residual_max = residual.lpNorm<Eigen::Infinity>();
    std::cout << "l2 residual is :" << residual_l2
              << " ,residual max norm is :" << residual_max << std::endl;
    Eigen::VectorXd error = Truesolution - possionproblem.getSolution();
    double error_l2 = error.norm();
    double error_max = error.lpNorm<Eigen::Infinity>();
    std::cout << "Error l2 norm is :" << error_l2
              << ", error max norm is :" << error_max << std::endl;

    // if (std::atoi(argv[2]) == 0)
    // {
    //     FILE *file1;
    //     file1 = fopen("output.txt", "w"); //
    //     “w”写文件覆盖原内容，不想覆盖则用“a” if (NULL == file1)
    //         printf("error");
    //     fprintf(file1, "%f\n%f\n%f\n%f\n", domain->getVertex(0)[0],
    //     domain->getVertex(1)[0], domain->getVertex(0)[1],
    //     domain->getVertex(3)[1]); fprintf(file1, "%d\n%d\n%d\n",
    //     m->getSegment()[0], m->getSegment()[1], 1); for (int i = 0; i <
    //     possionproblem.getSolution().size(); i++)
    //         fprintf(file1, "%.12f\n", possionproblem.getSolution()[i]); //
    //         循环输出数据并写入
    //     fclose(file1);                                                  //
    //     最后一定记得关闭文件 std::cout << "Data has been output" <<
    //     std::endl;
    // }
    // else if (std::atoi(argv[2]) == 1)
    // {
    //     FILE *file2;
    //     file2 = fopen("err.csv", "a");
    //     if (NULL == file2)
    //         printf("error");
    //     fprintf(file2, "%d %d %.12f\n", std::atoi(argv[1]),
    //     std::atoi(argv[2]), error_l2); std::cout << "Data has been output" <<
    //     std::endl;
    // }
    // else
    //     std::cout << "No output" << endl;
}
